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ABSTRACT 


There is considerable interest in analysing the effect of sample design on time series 
decomposition into components such as trend, seasonal and residual. One of the main 
deficiencies in past analysis has been a lack of a realistic model of an officially used time 
series decomposition method. This paper demonstrates a realistic approximation for X11 
and shows how such a model can be used to analyse the effect of sample design on time 
series decomposition. 


TIME SERIES DECOMPOSITION USING X11 AND SAMPLING ERRORS 
1. Introduction. 


There is considerable interest in the effects of sampling error on the components of time 
series decomposition, such as, the seasonally adjusted, trend and residual data. For 
example it is useful to have standard errors on the seasonally adjusted data, trend, seasonal 
and residual which represent the effects of the sampling process. An important question is 
what proportion of the observed variability in a time series is due to sampling error. Then 
there are design considerations, for example given a certain sample design and 
characteristics what is an "optimal" time series decomposition. One of the main 
deficiencies in past analysis has been providing a realistic model of an officially used time 
series decomposition method. This paper demonstrates a realistic model for X11(Shiskin 
et al, 1967) and shows how such a model can be used to analyse the effect of sampling 
design on time series decomposition. 


To give a appreciation of the complexity of the problem consider the simple time series 
model 


it y=t+s+e 
where 
1.2 t is the trend 
s is the seasonal 
é is the residual/irregular 


It is reasonable to assume the residual contains the following components 


1.3 @ =erwtes tens 
where 
1.4 erw is the "real world" noise 


es is the sampling error 
ens is the non-sampling error 


However the following sources of error may apply to the "true" decomposition : 


(i) The sampling error will affect the estimation of the trend, seasonal and "real 
world" noise. The sampling error may also contain "trend like" and 
seasonal behaviour. 
(il) The non-sampling error will affect the estimation of the trend, seasonal and 
"real world" noise. In addition the non-sampling error may well 
contain trend and seasonal components itself. 
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(iii) | Revision to the data will effect components to varying degrees. 


(iv) The components at the end of the data will be revised as more data becomes 
available due to the X11 method. That is assuming the central filters 
in X11 give the "true" components the estimates at the end of the 
data will only be an approximation. 


(v) It is assumed that X11 methodology can adequately estimate the trend, 
seasonal and residual. 


(vi) It is assumed that other components such as trading-day, moving holidays, 
and abrupt changes in the level or seasonal pattern are adequately 
estimated when they exist in the data. 


2. Model to analyse the effect of sampling error on time series decomposition 


The Australian Bureau of Statistics (ABS) has used the X11 method for all officially 
seasonally adjusted time series since 1967. Generally time series sample surveys have 
error estimates that are correlated. Hence the two main tools to enable an analysis of the 
effect of sampling errors on a time series are a realistic model of X11 and a covariance 
matrix associated with the sample design. 


3. A realistic model of X11. 


3.1 Previous models 

Young (1963) gave a basic linear approximation to X11. In general though this, and other 
subsequent work has been inadequate to accurately represent X11 as used in practice. 

The main area's where the previous models have been inadequate are :- 


(i) Studies have mainly concentrated on the central filters used in X11. These 
are not representative of the actual filters used by X11 (for example 
compare graph A3 with A5). 


(11) In the main the standard options available in X11 have been used. All other 
options have been ignored including the fact that X11 uses different 
filters for short time series (e.g. 5 years long). 


(iii) tis not recognised that logging the data and adjusting additively is not 
equivalent to the multiplicative option in X11 especially in terms of the 
arithmetic levels. 


(iv) | The modification for outliers used in X11 is usually completely ignored. 


(v) Effects such as trading-day, moving holidays, trend breaks, seasonal breaks 
etc. have been ignored. 


For example the often quoted paper by Cleveland et al. (1976), is clearly deficient on all of 
the above concerns. Other papers that are deficient in at least some of the considerations 
above include Hausman and Watson (1985), Maravall (1985), and Butter and Mourik 
(1990). A much more comprehensive treatment is given in the recently published working 
paper by Dagum, Chhab and Chiu (1993). 


3.2. Other Methods 

Some authors have assumed other theoretical models for the decomposition. These include 
SABL,STL (Cleveland et al.) which use moving medians and local regression 
respectively, BSM (Maravall et al. 1985) a basic structural model , and STM(Butter and 
Mourik, 1990) a structural model using a State Space representation and estimated using 
the Kalman filter. This approach ignores the fact that many major statistical organisations 
are using X11. In addition these alternative models are still deficient with respect to some 
of the points listed above. 


Professor James Durbin said in a speech in 1986 "...and a great deal of progress has been 
made in the theory of time series analysis. It would, therefore, have been natural to expect 
that improvements in methods of seasonal adjustment would have taken place at a rapid 
rate, and that the techniques in daily use today would have been revolutionised compared 
to the X11 method of twenty year ago. However this has not happened." Part of the 
problem even today is that while many new structural models have been proposed and 
offer promising alternatives to X11 they are not developed or tested to a form suitable for 
a large statistical organisation to decompose hundreds, or even thousands of time series 
with differing characteristics. 


3.3. Non-iterative version of X11. 


One of the design features of X11 is to provide an ad-hoc iterative procedure that 
converges to a relatively stable estimate quickly. A non-iterative procedure that gives very 
similar results to X11 is outlined below. A feature of this implementation is that a general 
matrix language (PROC IML in SAS) has been used. This has enabled compact code and 
all of the options available in X11 (and some that are not) to be included. 


3.4. Notation 
3.4.1 nby 1 vectors where n is the length of the data 


weights(0 no modification, | fully modified) 
seasonally adjusted data 


y original data 

i trend 

s seasonal 

@ residual/irregular 
w 

a 


3.4.2 nby n matrices : 
T A matrix to derive the trend from seasonally adjusted data. For example for 
a 7 term Henderson the matrix would look like: 


0.535 0.383 0.116 —0.034 0 ‘ 0 
0.289 0.410 0.294 0.061 -0.054 0 : 0 
0.034 0.275 0.399 0.287 0.058 —0.053 0 0 
—0.059 0.059 0.294 0.412 0.294 0.059 -0.059 0 
T= 0 
0 
0 
0 : : : : 
0 —0.034 0.116 0.383 0.535 
3.4.3 SI A matrix to compute the seasonal factors from the detrended data. 
For example a 3 by 5 moving average. 
3.4.4 $2 A matrix to correct for levels.. For example a 2 by 12 moving average. 
Let 
3.4.5 S=(—S$2)*S1 


Then for a simple additive model 
3.4.6 y=tt+st+e 


we have three equations 


3.4.7 t =T*(y-35) 
5 =S*(y-t) 
@=y-t-s 

Hence 

3.4.8 t=T*G-S*(-f)) 


(I-T*S)*t =(T-T*S)*y 
Providing that the inverse exists we have 
3.4.9 t =(-T*S)7 *(T-T*S)*y 


Hence the component's ¢ ,5,@ and the adjusted data @ can be directly computed from J. 


3.4.10 i=T*y where T =(1-T*S)1*(T-T*S) 
F=S*7P where S=S-S*T 
e=E*y where E=I-T-S 
G=A*y where A=(-S) 


3.5. Multiplicative adjustment 
If 


3.5.1 yor *s ¥e 
then logging 3.3.1 gives an additive model in logs 
3.5.2 log(V) = log(t )+log(S) + log(é) 


Hence estimates of the components are given by 


3.5.3 7 =exp(T * log(¥)) 
s =exp(S * log(y)) 
é =exp(E * log(y)) 


Unfortunately this leads to geometric moving averages being applied. This means that 
levels are maintained in geometric rather than arithmetic terms (as say multiplicative X11 
does). It is interesting to note that the symmetric Henderson weights do have the property 
that the weights applied to logged data and then un-logged are very close to the weights 
applied directly to the un-logged data (because Henderson moving averages leave 
polynomial trends of up to cubics unchanged). The end weights and the moving averages 
used to obtain the seasonal factors from the de-trended data do not have this property. 
Hence, if a seasonal series is adjusted using multiplicative X11 and compared to the same 
adjustment logging the data, applying an additive adjustment and then un-logging the data 
a systematic bias is observed in say the trend levels. There are many ways to attempt to 
correct for this bias. A simple and effective way is to correct the seasonal factors obtained 
with the log additive model for multiplicative level bias. 


That is 


3.5.4 sbia¥ = exp(S * log(¥)) 
S = sbias (S2 * sbias )) 


3.6. Modifications for extremes 

An important part of X11 is to deal with outliers. To accurately represent X11 (and to 
provide good estimates of the components) outliers must be modified. A simple and 
effective method is to use the concept of a weight applied to the data. 


That is, let w be a vector of weights to apply, where w(i) =0 gives no modification and 
w(i) = 1 gives full modification, and 0 < w(i) < 1 gives partial modification. Then it can be 
shown that the matrices to derive the trend, seasonal and residual can be modified to allow 
for outliers as shown below. 


Assuming T,S and E have been computed as outlined above then they can be modified to 
take account of outliers as follows 


3.6.1 E=(I-(I-E) * diag(w))1* E 
T=T (I—diag(W) * B) 
S=S * (I—diag(W) * E) 


The weights w can be taken from X11 or directly computed using a method similar to 
X11. In theory the computation of the weights would need to be taken into account for 
standard errors on the estimates. This has not been attempted in this paper. 


Several alternative methods of handling outliers are currently being investigated. These 
include modifications of the form 


3.6.2 y-r 

where A is of the form 

203 R=g*(+g*E!* EB)! * EF! * E*F where 0<g<oo 
or with suitable G 

3.6.4 R=g*([+9*G! * BE! *E*G)! *G! * El * E*F 


One interesting application is the choice of G to enable outliers to receive more smoothing 
rather that being removed or reduced in the data. 


3.7. Trading-day, moving holidays and other influences 

Others types of components such as trading-day, moving holidays and abrupt changes in 
the level or seasonal pattern can be incorporated into the methodology outlined above. 
For example letting X be an n by k matrix where k is the number of parameters to be 


estimated and X appropriately formulated for trading-day (Young 1965) or moving 
holidays. Then these components can be incorporated as 


S74 D=(1-H*(I-E))"' *H*E 
where H=X * (X/ * X)! # X’, 


The other components are modified as 


bere) Paper D) 
S=S*(I- D) 
E=f2TH8=D 


3.8. Extensions to X11 

The basic X11 algorithm can easily be extended to include other models, for example it is 
a simple extension to allow for regression estimates of the trend and seasonal components. 
The residual component could be allowed to follow a moving average process 
(Box-Jenkins). For example for the basic model 


3.8.1 y=t+s+e 


the residual @ could follow a moving average process of order 1 


3:8.2 e(t) = e€(t) -@ * e(t— 1) where €~N(0, 07) 


letting 
1000 0. 
6100 0 
3.8.3 Galion we oe 0 
0 
69 1 0 
0 00-61 
we have 
3.8.4 € =(I-E*(I-G))"| * E*F 


Hence 6 could be estimated by least squares by minimising €’ * € or with appropriate 
modifications by maximum likelihood. 


3.9. Implementation 

The vector and matrix formulation given above is ideally suited to Proc IML in SAS, and 
other matrix oriented packages. All computations have been done using IML. The 
computations for the Henderson moving averages have been done using algorithms for the 
central and surrogate weights (used at the ends of the data) to make the method more 
general. 


As an example Australian Monthly Retail Turnover has been used. The data from the 
1992 annual re-analysis with the same prior correction factors has been used. The span of 
the analysis has been restricted to 7/83 to 6/92 due to the large amount of computer 
memory required. 


The seasonal factors are computed using the above methodology from 7/83 to 6/92 and the 
forward factors computed using the same method as in X11. The adjusted series is 
compared to the published figures (April 1993) in graph Al. The movements are 
compared in Graph A2. For this example the differences are minor. Several other time 
series have been looked at with similar results. More comprehensive testing is envisaged 
in the future. 


Graphs A3 to Al4 shows the weights for some selected point in time for the trend, 
seasonal and irregular components. Some of the features are that the weighting patterns at 
the ends of the data are completely different to those at the centre (e.g. compare A3 with 
A5). The modified weighting patterns may be similar to the unmodified in some cases, 
and different in others, depending on the weights and the location of the extremes (e.g. 
compare A7 and A8). 


Graph A9 shows T graphically for the case of a 9 term Henderson, 3x5 seasonal and 
modifications for extremes for the Monthly Retail Turnover data. The left hand axis is the 


actual weights, the bottom axis is the time point at which the filter given by the right axis 
is applied. 


While not shown in this paper the spectral properties of the filters outlined above can be 
extensively analysed using the gain and phase outputs of a linear filter. 


4. A covariance matrix associated with the sample design 


The computation of the covariance matrix associated with the sampling design is mainly a 
sampling problem. There seems to be three approaches possible. These are : 


(i) Compute the whole covariance matrix directly using the sampling design. 
There would usually be computational difficulties in computing such 
a matrix going back several years. 


(ii) Find a model to estimate the covariances. This was the approach taken by 
Steel and De Mel (1987) for the Australian Labour Force Data where a 
"geometric decay" model is used. 


(iii) | Make up the covariance matrix using an educated guess for the model. 


It should be noted that any covariance matrix used must be positive semi-definite to ensure 
zero Or positive variances. 


The analysis would be considerably complicated if the sampling error, non-sampling error 
and "real world" error are correlated. There is also the somewhat philosophical problem of 
looking at the "true" observed values as fixed, versus their own distribution, that is are the 
sample totals fixed or random variables. Given these tools the type of analysis that can be 
achieved is outlined below. 


The basic theoretical framework is outlined in section 3. Ignoring a correction for levels, 
the components of a multiplicative time series analysis, with modification for extremes 
(using the modifications in 3.4), is given by 


4.1 7 =exp(T * log) 
5 =exp(S * log(y)) 
é =exp(E * log(y)) 


5. Additive verses multiplicative standard errors. 
The published standard errors are usually given as an additive standard error. Given the 
multiplicative option in X11 is almost always used this means that the basic model is given 


by 


5.1 y=t *5 *erw* ens tes 


This presents some problems, namely: 
(i) Logging the data does not give a model which is easy to analyse. 


(il) Such a model cannot be totally reasonable since it implies negative values 
are possible in data that cannot be usually negative. 


(iii) The standard errors are usually relative to the level of the data. For example 
the published levels, standard errors and percentage standard errors 


for Australian Monthly Retail Turnover is given in the table 

below 

Sid TABLE 1 

Date Level Standard Error % Standard Error 

02/92 7,106.2 78.9 1.11 
03/92 7,475.9 73.6 0.98 
04/92 7,694.5 fie) 0.96 
05/92 tal 2 74.8 0.96 
06/92 7,547.9 63.4 0.84 
07/92 7,819.7 62.4 0.8 
08/92 7,461 59.4 0.8 
09/92 7,745.6 60.4 0.78 
10/92 8,252.3 64.4 0.78 
11/92 8,126.9 64.8 0.8 
12/92 10,637.7 Ol 0.86 
01/93 7,789.5 67.2 0.86 
02/93 7,108 995 0.84 
03/93 7,831.6 66.2 0.85 
04/93 7,905.4 66.5 0.84 


Because of the above problems it is more reasonable to have the covariance matrix in 
terms of a multiplicative error. That is 


5.3 y=t *5 *erw*ens *(1+65) 
6. Variance on the estimated components 
Assuming a relative covariance matrix C, then variances of the components are given by 
i x vs 
6.1 02,(t ) =diag(T * C*T ) 
/ 


62,(5) =diag(S*C*S ) 
10 


02,(é) =diag(E* C*E ) 


Hence 95 per cent confidence bounds (assuming +2 sigma is 95%) for the components are 
given by 


6.2 t *(1+2*0,5(f)) 
S *(14£2 * 0,s(5)) 
@ *(142*0,,(e)) 


From these formulas additive errors can be numerically computed. It should be noted that 
these estimates are highly correlated, and hence cannot be used to give simultaneous 
confidence bounds on several time points. In addition these formulas will only provide 
estimates of the uncertainly due to the sampling error. 


To give an example Monthly Retail Turnover data has been analysed. It should be noted 
that this analysis is for demonstration purposes only, and is not necessarily realistic. In 
practice the covariance matrix associated with the sample design would include explicit 
allowance for the rotation used, sample redesign and other known sampling design 
characteristics. The sample design area of the ABS is currently researching such 
covariance matrices for the ABS sample surveys. 


An examination of table 1 above and the standard errors on the movements given in the 
publication a standard error of about | per cent and a high correlation between successive 
time points might be a reasonable model for the covariance matrix. A possible model 
might be an AR(1) model of the form 


6.3 es; =P *es;-1+€; where €, iid with variance 67 


it can be shown that this has the covariance matrix 


1 ) 7 . po”! 
p 1p pie 
6.4 o 
I-p : , ; : 
Pie g epe ae Gp 
prt ; 7 te) 1 


Hence letting p = 0.8 and o? = 0.0001 * (1 — 0.87) gives a 1 per cent standard error and 
correlation at lag 1 of 0.8 (an alternative model might have been a MA(1)). 


Graphs 1,2 ,3 and 4 show the standard error of the trend, seasonal factors, irregular and 
seasonally adjusted components as a proportion of the original standard error (1%). For 
example Graph | shows that the trend from X11 is about 104 per cent of the standard error 
on the original data at the end of the data, and about 90 per cent in the middle. These 
results would change if the options in X11 are varied, the covariance matrix is different 
(either different parameter values, or different model). It should be noted that the trend 
standard error is for the X11 trend, which in the case of Retail is different from the 
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published trend. This is because ABS currently uses a 9 term Henderson for the seasonal 
adjustment, while the published trend is uses a 13 term Henderson, and no allowance is 
made for outliers in the published trend. 


It is relatively straightforward to modify the procedure to allow for these differences. The 
approximate standard error of the published trend as a proportion of the original standard 
error is shown in graph 5. If the proportions of the trend in graph | and graph 5 are 
compared it will be noticed that it is lower in graph 5 (NB 13 vs 9 Henderson). However 
this does not however imply that estimation of the trend without modifications for 
extremes is superior. There are two competing factors in estimating the trend (or any other 
component) "bias" and "variance". That is the trend produced without modifications for 
extremes may have a lower variance but be also very biased in producing what is deemed 
a reasonable trend. 


A similar situation is faced at the end of the data. It has been noticed for some time series 
that the trend produced at the end of the data seems to be biased when compared to the 
final trend. Such bias could be eliminated by applying the same criteria that are used to 
compute the central Henderson weights at the end of the data. The result would be a much 
higher variance on the trend at the end of the data. 


Comparing these results with the work of Steel and De Mel (op. at.) shows that while there 
are broad similarities in the magnitude, it is clear that there is considerably more 
complexity in the filters actually used in X11 to obtain the components. Some interesting 
features are the pronounced rises in the early years of the proportion for the seasonal 
factors. This is due to a predominance of full modification for outliers in those months 
(compare with graph 6). 


The variability of the residual due to the sampling process can be compared with the 
estimate of variance for the observed residual. This includes "real world" error, sampling 
error and non-sampling error, and (assuming a constant variance), is given by : 


6.5 02 =F Li, (e,-1)? 


Hence the percentage of volatility of the levels of Monthly Retail Turnover due to the 
sampling process is approximated by 


6.6 100 « 2@ 


oz(@) 


Graph 7 shows this percentage over time. In practice, the "real world" variance is almost 
certainly changing over time and there is plenty of empirical evidence that it is seasonal. 
It is a moot point whether the proportion should be to the actual observed residual or the 
residual modified for outliers. In the latter case the proportion would be much higher. 
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Note: 9 term Henderson, 3x5 seasonal, modification for extremes. 
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Note: 9 term Henderson, 3x5 seasonal, modification for extremes. 


MONTHLY RETAIL TURNOVER 
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Note: 9 term Henderson, 3x5 seasonal, no modification for extremes. 
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MONTHLY RETAIL TURNOVER 
Variance of sampling verses observed error 


GRAPH 7 
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Note: 9 term Henderson, 3x5 seasonal, modification for extremes. 
Observed variance assumed constant. 


7. Variance of the movements. 


These are estimated in a similar fashion to that of the levels shown in section 4. In this 
case define 


7A DpDi= T-—BT where B, the shift operator, shifts the rows of a matrix 
down one. 

D’ =S-BS 

Dé =E-BE 


Then the approximate variances on the movements due to the sampling process is given by 


12 o2,((t — Bt )/Bt ) = diag(D' * C * D") 
62,((§ — BS)/BS) = diag(D’ * C* D") 
o2,((é — Bé)/Be) = diag(D* * C * D”) 


Again, the variance of the change in the irregular component can be compared with the 
estimated variance of the observed change in the residual. This is given by 


23 02((é —Bé)/Be) = + YL (<5- 1)? 


Hence the percentage of volatility of the movements of Monthly Retail Turnover due to 
the sampling process is approximately given by 


Ges((€-BE/B2) 
02 ((€-Bé)/Bé) 


7.4 100 * 


Graph 8 shows this percentage over time. 


8. Can better estimates of the "real'' data be computed? 


There are several approaches that can be used to provide better estimates of the "real" data. 
Unfortunately because there are so many sources of error any method proposed can 
always be criticised because of the assumptions used. 


One approach is to use spectral analysis. For example if it could be ascertained, either 
with prior knowledge or empirical investigation, that the contribution of the variance due 
to sampling error was concentrated in a certain bandwidth and the "real world" and 
non-sampling error in a different bandwidth, then an appropriate filter could be designed 
using spectral analysis to remove/dampen the sampling error and leave the other error 
unchanged. More formally if the model for the observed data is 


8.1 yor +s +es +er 
where 
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8.2 es is the sampling error and 
er is the rest of the error. 


Then if one could find a filtering matrix F such that 
8.3 er =F * (es +er)+€ 


then it can be shown that a better estimate of the original data is given by 
Aa Aa -l Aa 
8.4 i=(1-(E+F*u-B) sPxu-n) 3 


8.5 =R*y 


There are several ways a suitable F matrix could be estimated. For example, if the 
sampling error is assumed to follow an autoregressive model then it can be written as 


8.6 Ges =€ 
Then if / depended on parameters they could be estimated by minimising 


8.7 ¥ *(G*R) *(G*R)*F 


9. Sample design using time series characteristics. 


If it is required to have an "optimal" sample allocation with respect to the decomposed 
data then clearly time series characteristics of the data must be taken into account. For 
example Australian Monthly Retail Turnover has allocation goals of minimising standard 
errors on movement and level for the original data. Generally if estimating the change in 
level between two time periods estimates of maximum precision are obtained by retaining 
the same sample on both occasions. For average or total values over a number of surveys 
non-overlapping units are selected. Clearly the trend level and movements are a weighted 
average of several survey values. If the goal was for a optimal sample for the trend 
component from the decomposition then the sample design/allocation may well be 
different. 


In the example for Australian Monthly Retail Turnover the observed variance was 
assumed to be constant. In practice the variance may well be changing over time and is 
often observed to be seasonal. 


It is possible to compute a time varying variance in a similar manner to that of the levels. 
For example, using a 5 year moving average and constant seasonality a model for the 
variance might be ( where V is derived in a similar way to that of the levels in section 3) 


9.1 6°(@)=V¥*e- 
where 
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In which case 9.1, could be used to assist in deriving an optimal sample allocation. 


It should also be pointed out that autocorrelation in the residual due to sample design error 
(if strong enough) could be used to determine optimal filters for the time series 
decomposition. 


10. Conclusions. 


This paper demonstrates that a realistic approximation of X11, the widely used time series 
decomposition package, is possible. This provides a tool to allow extensive analysis of 
sample design and characteristics on time series decomposition. In particular, given a 
model for the covariance matrix of the sample design, and making certain assumptions the 
paper shows that 

standard errors can be computed on components such as trend, seasonal and residual at all 
time points, and the proportion of variance due to sampling error can be estimated. 


This paper only outlines the tools required for the analysis and further work needs to be 
completed on practical applications for individual surveys. 


Currently, sampling design concentrates only on the original data. Given that the 
Australian Bureau of Statistics is giving more emphasis on the "trend" estimates it is not at 
all obvious that an "optimal" sample design for the original data will be "optimal" for the 
trend as estimated by the Australian Bureau of Statistics. This is clearly an area where 
further work is required to integrate the sample design and decomposition process. 
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